# Import PyDicom before using PyDicom functions
# DICOM (Digital Imaging and Communications in Medicine) is the standard protocol for
# the management and transmission of medical images and related data, and is used by
# many healthcare facilities
# PyDicom is a pure Python package for working with DICOM files
import pydicom
# The `pathlib` module is similar to the `os.path` module, but `pathlib` provides a more
# advanced and convenient interface than `os.path`
# It is possible to use `pathlib` to represent file paths as specialized `Path` objects
# instead of plain strings
from pathlib import Path
# SimpleITK is a simplified programming interface to the algorithms and data structures
# of the Insight Toolkit (ITK) for segmentation, registration and advanced image analysis
# SimpleITK supports bindings for multiple programming languages including C++, Python, R,
# Java, C#, Lua, Ruby and TCL
# Combining SimpleITK’s Python bindings with the Jupyter notebook web application creates
# an environment which facilitates collaborative development of biomedical image analysis
# workflows
# In this project, SimpleITK provides the ability to automatically detect and read all DICOM
# files so that the user does not have to manage file reading or slice sorting
import SimpleITK as sitk
import numpy as np
import torch
from torchvision.utils import make_grid
import matplotlib as mpl
import matplotlib.pyplot as plt
from datetime import datetime
from functools import wraps
import itertools
import math
import random
import reprlib
import sys
%matplotlib inline
XINHUI = "#7a7374"
XUEBAI = "#fffef9"
YINBAI = "#f1f0ed"
YINHUI = "#918072"
figure_size = (16, 9)
custom_params = {
"axes.axisbelow": True,
"axes.edgecolor": YINBAI,
"axes.facecolor": XUEBAI,
"axes.grid": True,
"axes.labelcolor": XINHUI,
"axes.spines.right": False,
"axes.spines.top": False,
"axes.titlecolor": XINHUI,
"figure.edgecolor": YINBAI,
"figure.facecolor": XUEBAI,
"grid.alpha": 0.8,
"grid.color": YINBAI,
"grid.linestyle": "--",
"grid.linewidth": 1.2,
"legend.edgecolor": YINHUI,
"patch.edgecolor": XUEBAI,
"patch.force_edgecolor": True,
"text.color": XINHUI,
"xtick.color": YINHUI,
"ytick.color": YINHUI,
}
mpl.rcParams.update(custom_params)
reprlib_rules = reprlib.Repr()
reprlib_rules.maxother = 250
sys.path.append("../")
from Modules import *
# The dataset used in this practice project is a very small subset of CT images extracted from
# the Cancer Imaging Archive (TCIA), which contains intermediate slices of all CT images
# in which valid age, mode, and contrast markers can be found
# The dataset is provided by the Kaggel Datasets called CT Medical Imaging
# (https://www.kaggle.com/datasets/kmader/siim-medical-images),
# license type is Attribution 3.0 Unported (CC BY 3.0)
# https://creativecommons.org/licenses/by/3.0
# The link to the TCIA archive of the full dataset is
# https://wiki.cancerimagingarchive.net/display/Public/TCGA-LUAD
dir_path = "../Datasets/Kaggle - CT Medical Images/dicom_dir/"
sample_dcm = "ID_0000_AGE_0060_CONTRAST_1_CT.dcm"
# The main function of PyDicom to read and parse DICOM files is `read_file`
dicom_file = pydicom.read_file(dir_path + sample_dcm)
tabulation = Form_Generator()
tabulation.heading_printer("Initial understanding of DICOM file")
statements = [
"""
dir_path = "../Datasets/Kaggle - CT Medical Images/dicom_dir/"
sample_dcm = "ID_0000_AGE_0060_CONTRAST_1_CT.dcm"
dicom_file = pydicom.read_file(dir_path + sample_dcm)
"""
]
tabulation.statement_generator(statements)
variables = ["dicom_file"]
values = [str(reprlib_rules.repr(dicom_file))]
tabulation.variable_generator(variables, values)
expressions = [
"dicom_file[0x0028, 0x0010]",
"dicom_file[0x0028, 0x0011]",
"dicom_file[0x0018, 0x0015]",
"dicom_file.Rows",
"dicom_file.Columns",
"dicom_file.BodyPartExamined",
"dicom_file.keys()",
"dicom_file.values()",
"dicom_file.dir()",
'dicom_file.dir("Image")',
]
results = [
str(dicom_file[0x0028, 0x0010]),
str(dicom_file[0x0028, 0x0011]),
str(dicom_file[0x0018, 0x0015]),
str(dicom_file.Rows),
str(dicom_file.Columns),
str(dicom_file.BodyPartExamined),
str(reprlib_rules.repr(dicom_file.keys())),
str(reprlib_rules.repr(dicom_file.values())),
str(reprlib_rules.repr(dicom_file.dir())),
str(reprlib_rules.repr(dicom_file.dir("Image"))),
]
tabulation.expression_generator(expressions, results, 1)
Initial understanding of DICOM file +-------------------------------------------------------+ | Statement | +-------------------------------------------------------+ | dir_path = "../Datasets/Kaggle - CT Medical | | Images/dicom_dir/" | | sample_dcm = "ID_0000_AGE_0060_CONTRAST_1_CT.dcm" | | dicom_file = pydicom.read_file(dir_path + sample_dcm) | +-------------------------------------------------------+ +------------+------------------------------------------------+ | Variable | Value | +------------+------------------------------------------------+ | dicom_file | Dataset.file_meta | | | ------------------------------- | | | (0002, 0000) File Meta Information Group | | | Length UL: 194 | | | (0002, 0001) Fil...Group Length | | | UL: 524296 | | | (7fe0, 0010) Pixel Data | | | OW: Array of 524288 elements | +------------+------------------------------------------------+ +-----------------------------+-------------------------------+ | Expression | Result | +-----------------------------+-------------------------------+ | dicom_file[0x0028, 0x0010] | (0028, 0010) Rows | | | US: 512 | | dicom_file[0x0028, 0x0011] | (0028, 0011) Columns | | | US: 512 | | dicom_file[0x0018, 0x0015] | (0018, 0015) Body Part | | | Examined | | | CS: 'CHEST' | | dicom_file.Rows | 512 | | dicom_file.Columns | 512 | | dicom_file.BodyPartExamined | CHEST | | dicom_file.keys() | dict_keys([(0008, 0000), | | | (0008, 0005), (0008, 0008), | | | (0008, 0016), (0008, 0018), | | | (0008, 0020), (0008, 0021), | | | (0008, 0022), ...095, 0010), | | | (0097, 0000), (0097, 0010), | | | (0099, 0000), (0099, 0010), | | | (7003, 0000), (7003, 0010), | | | (7fe0, 0000), (7fe0, 0010)]) | | dicom_file.values() | dict_values([(0008, 0000) | | | Group Length | | | UL: 430, (0008, 0005) | | | Specific Character Set | | | CS:...up Length | | | UL: 524296, | | | (7fe0, 0010) Pixel Data | | | OW: | | | Array of 524288 elements]) | | dicom_file.dir() | ['AccessionNumber', | | | 'AcquisitionDate', | | | 'AcquisitionNumber', | | | 'AcquisitionTime', | | | 'BitsAllocated', | | | 'BitsStored', ...] | | dicom_file.dir("Image") | ['ImageDimensions', | | | 'ImageFormat', | | | 'ImageGeometryType', | | | 'ImageLocation', | | | 'ImageOrientation', | | | 'ImageOrientationPatient', | | | ...] | +-----------------------------+-------------------------------+
# DICOM data has an attribute called `pixel_array` that provides more useful pixel data
# for uncompressed images that can be passed to the graphics library for viewing
# To use this attribute, the system must have the NumPy numeric package installed,
# since `pixel_array` returns a NumPy array
ct = dicom_file.pixel_array
tabulation = Form_Generator()
tabulation.heading_printer("Getting pixel data from DICOM file")
statements = ["ct = dicom_file.pixel_array"]
tabulation.statement_generator(statements)
variables = ["ct"]
values = [str(reprlib_rules.repr(ct))]
tabulation.variable_generator(variables, values)
expressions = ["ct.shape"]
results = [str(ct.shape)]
tabulation.expression_generator(expressions, results)
Getting pixel data from DICOM file +-----------------------------+ | Statement | +-----------------------------+ | ct = dicom_file.pixel_array | +-----------------------------+ +----------+------------------------------------------------+ | Variable | Value | +----------+------------------------------------------------+ | ct | array([[0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | ..., | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0]], dtype=uint16) | +----------+------------------------------------------------+ +------------+------------+ | Expression | Result | +------------+------------+ | ct.shape | (512, 512) | +------------+------------+
def image_display(image, ax, title, cmap):
ax.imshow(image, cmap)
ax.grid(False)
ax.set_title(title, loc="center", pad=10)
x_ticks = list(range(0, image.shape[1], 50))
y_ticks = list(range(0, image.shape[0], 50))
ax.set(xticks=x_ticks, xticklabels=x_ticks,
yticks=y_ticks, yticklabels=y_ticks)
ax.set_xlim(left=0)
ax.set_ylim(top=0)
ax.minorticks_on()
return ax
# Path classes are divided between pure paths, which provide purely computational
# operations without I/O, and concrete paths, which inherit from pure paths but also
# provide I/O operations
path_object = Path(dir_path)
# `PurePath.name` returns a string representing the final path component, excluding
# the drive and root directory (if any)
# When a path points to a directory, `Path.iterdir` generates a path object of the directory's
# contents
# `Path.is_file` returns True if the path points to a normal file or to a symbolic link
# to a normal file, False if it points to another type of file
random_dicom_path = random.choice(
[
file.name
for file in path_object.iterdir()
if file.is_file() & (pydicom.read_file(file).BodyPartExamined != "CHEST")
]
)
random_dicom_file = pydicom.read_file(dir_path + random_dicom_path)
plt.rcParams["figure.figsize"] = (
figure_size[0] / 3 * 2, figure_size[1] / 4 * 5)
fig, axs = plt.subplots(nrows=2, ncols=2)
image_display(
ct,
axs[0, 0],
f"{dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'gray'",
cmap="gray",
)
image_display(
ct,
axs[0, 1],
f"{dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'bone'",
cmap="bone",
)
image_display(
random_dicom_file.pixel_array,
axs[1, 0],
f"{random_dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'gray'",
cmap="gray",
)
image_display(
random_dicom_file.pixel_array,
axs[1, 1],
f"{random_dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'bone'",
cmap="bone",
)
fig.suptitle(
"Visual Comparison of CT Medical Image Display Using Different Colormaps",
fontsize="x-large",
x=0.5,
y=0,
)
plt.tight_layout()
plt.show()
def annotation_color(cmap):
cmap_color = mpl.colormaps[cmap]
return cmap_color(0.85)
def first_max(list):
max_element = list[0]
for element in list[1:]:
if len(element) > len(max_element):
max_element = element
return max_element
def text_aligner(text):
line_list = text.split("\n")
max_length = len(first_max(line_list))
line_list = [
(":" + " " * (max_length - len(line))).join(line.split(":"))
if len(line) < max_length
else line
for line in line_list
]
text = "\n".join(line_list)
return text
def plane_information_annotator(func):
@wraps(func)
def wrapper(ax, cmap, dicom_file, fontsize_adjustment):
plane = func(ax, cmap, dicom_file, fontsize_adjustment)
text = f"Image: {dicom_file.Columns} x {dicom_file.Rows}"
text += f"\nImage Plane: {plane}"
text += f"\nPatient Position: {dicom_file.PatientPosition}"
text += f"\nSlice Thickness: {dicom_file.SliceThickness:.1f} mm"
text += f"\nSlice Location: {dicom_file.SliceLocation:.1f} mm"
ax.text(
0.575,
0.825,
text_aligner(text),
transform=ax.transAxes,
fontsize=9 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
return wrapper
def orientation_annotator(func):
@wraps(func)
def wrapper(ax, cmap, dicom_file, fontsize_adjustment):
plane, top, right, bottom, left = func(dicom_file)
ax.text(
0.4875,
0.95,
top,
transform=ax.transAxes,
fontweight="heavy",
fontsize=11 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
ax.text(
0.4875,
0.025,
bottom,
transform=ax.transAxes,
fontweight="heavy",
fontsize=11 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
ax.text(
0.95,
0.4875,
right,
transform=ax.transAxes,
fontweight="heavy",
fontsize=11 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
ax.text(
0.025,
0.4875,
left,
transform=ax.transAxes,
fontweight="heavy",
fontsize=11 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
return plane
return wrapper
def image_orientation(func):
@wraps(func)
def wrapper(*args, **kwargs):
Left, Right, Anterior, Posterior, Head, Feet = "L", "R", "A", "P", "H", "F"
axial_plane = [Anterior, Right, Posterior, Left]
coronal_plane = [Right, Head, Left, Feet]
sagittal_plane = [Anterior, Head, Posterior, Feet]
plane = func(*args, **kwargs)
# Patient Position (0018,5100) specifies the position of the patient relative to
# the space of the imaging device and is used for annotation purposes only,
# but does not provide a precise mathematical relationship between the patient and
# the imaging device
patient_position = dicom_file.PatientPosition
if plane == "Axial":
term_list = ["S", "DR", "P", "DL"]
top = 4 - term_list.index(patient_position[2:])
if patient_position.startswith("HF"):
right = top + 1
left = right + 2
else:
left = top + 1
right = left + 2
bottom = top + 2
top, right, bottom, left = [
axial_plane[index] if index < 4 else axial_plane[index - 4]
for index in [top, right, bottom, left]
]
elif plane == "Coronal":
term_list = ["DL", "DR"]
top = term_list.index(patient_position[2:]) * 2
if patient_position.startswith("AF"):
right = top + 1
left = right + 2
else:
left = top + 1
right = left + 2
bottom = top + 2
top, right, bottom, left = [
coronal_plane[index] if index < 4 else coronal_plane[index - 4]
for index in [top, right, bottom, left]
]
elif plane == "Sagittal":
term_list = ["LF", "RF"]
if patient_position.endswith("S"):
top = 0
right = term_list.index(patient_position[:2]) * 2 + 1
left = right + 2
else:
top = 2
left = term_list.index(patient_position[:2]) * 2 + 1
right = left + 2
bottom = top + 2
top, right, bottom, left = [
sagittal_plane[index] if index < 4 else sagittal_plane[index - 4]
for index in [top, right, bottom, left]
]
else:
top, right, bottom, left = np.repeat("", 4)
return plane, top, right, bottom, left
return wrapper
@plane_information_annotator
@orientation_annotator
@image_orientation
def get_plane_information(dicom_file):
# The DICOM standard defines a patient-oriented reference coordinate system (RCS) that
# enables the user to measure the position and orientation of the image relative to
# the patient
# Image Orientation (Patient) (0020,0037) specifies the cosine of the orientation of
# the first row and column relative to the patient, which should be provided as a pair,
# where the row values for the x, y, and z axes are followed by the column values for
# the x, y, and z axes
# The orientation of the axes is determined solely by the orientation of the patient, i.e.,
# the RCS and Image Orientation (Patient) (0020,0037) values specify the orientation of
# the image frame rows and columns
IOP = dicom_file.ImageOrientationPatient
row_xyz = IOP[:3]
column_xyz = IOP[3:]
# `numpy.cross` returns the cross product of two vector arrays
# The geometric interpretation of the cross product is that two vectors determine a plane,
# and the cross product points in a direction different from both vectors
# The Patient-Based Coordinate System is a right-handed system, i.e., the vector
# cross product of a unit vector along the positive x-axis and a unit vector along
# the positive y-axis is equal to a unit vector along the positive z-axis
plane = [abs(round(x)) for x in np.cross(row_xyz, column_xyz)]
# The Image Type (0008,0008) identifies important image identification features
# and it returns the Pixel Data Characteristics, Patient Examination Characteristics,
# Modality Specific Characteristics, and Implementation Specific Identifiers
# While the third value (Modality Specific Characteristics) should identify any
# specialization specific to the Image Information Object Definition (IOD), it
# and the values that follow it are not mandatory, unlike the first two values,
# and therefore some DICOM data does not have this value
if plane[0] == 1 and plane[1] == 0 and plane[2] == 0:
return "Sagittal"
elif plane[0] == 0 and plane[1] == 1 and plane[2] == 0:
return "Coronal"
elif plane[0] == 0 and plane[1] == 0 and plane[2] == 1:
return "Axial"
else:
"Unknown"
def get_patient_information(ax, cmap, dicom_file, fontsize_adjustment, patien_age):
text = f"Patient ID: {dicom_file.PatientID}"
text += f"\nPatient's Sex: {dicom_file.PatientSex}"
try:
patien_age = dicom_file.PatientAge.strip(str(0))[:-1]
except AttributeError:
patien_age = patien_age
text += f"\nPatient's Age: {patien_age} y"
text += f"\nModality: {dicom_file.Modality}"
text += f"\nBody Part Examined: {dicom_file.BodyPartExamined}"
ax.text(
0.04,
0.825,
text_aligner(text),
transform=ax.transAxes,
fontsize=9 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
def get_date_information(ax, cmap, dicom_file, fontsize_adjustment):
date = datetime.strptime(dicom_file.StudyDate,
"%Y%m%d").strftime("%Y.%m.%d")
text = f"Study Date: {date}"
ax.text(
0.04,
0.05,
text,
transform=ax.transAxes,
fontsize=7 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
def get_manufacturer_information(ax, cmap, dicom_file, fontsize_adjustment):
text = f"Manufacturer: {dicom_file.Manufacturer}"
ax.text(
0.525,
0.05,
text.rjust(36),
transform=ax.transAxes,
fontsize=7 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
def DICOM_2D_image(ax, cmap, dicom_file, title, fontsize_adjustment=0, patien_age="__"):
ax.imshow(dicom_file.pixel_array, cmap)
ax.grid(False)
ax.set_title(title, loc="center", pad=10,
fontsize=12 + fontsize_adjustment)
ax.set(xticks=[], yticks=[], frame_on=False)
get_patient_information(ax, cmap, dicom_file,
fontsize_adjustment, patien_age)
get_plane_information(ax, cmap, dicom_file, fontsize_adjustment)
get_date_information(ax, cmap, dicom_file, fontsize_adjustment)
get_manufacturer_information(ax, cmap, dicom_file, fontsize_adjustment)
fig, axs = plt.subplots(2, 2, figsize=(
figure_size[0] / 3 * 2, figure_size[1] / 4 * 5))
DICOM_2D_image(
axs[0, 0],
"gray",
dicom_file,
f"{dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'gray'",
)
DICOM_2D_image(
axs[0, 1],
"bone",
dicom_file,
f"{dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'bone'",
)
DICOM_2D_image(
axs[1, 0],
"gray",
random_dicom_file,
f"{random_dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'gray'",
)
DICOM_2D_image(
axs[1, 1],
"bone",
random_dicom_file,
f"{random_dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'bone'",
)
fig.suptitle(
"Visualization of CT Medical Image Display with Supplementary Information",
fontsize="x-large",
x=0.5,
y=0,
)
plt.tight_layout()
plt.show()
# The dataset used in this practical project is the MRI DICOM dataset of the head of
# a 52-year-old normal male mathematics professor
# The subject of the experiment is the author, who suffers from small vertical strabismus
# (hypertropia), which is a misalignment of the eyes, which can be seen in this dataset
# This dataset was provided by Lionheart, William R.B. in 2015, available online
# (https://zenodo.org/record/16956 or http://doi.org/10.5281/zenodo.16956),
# license type is Attribution-ShareAlike 4.0 International (CC BY-SA 4.0)
# https://creativecommons.org/licenses/by-sa/4.0
path_to_head_mri = Path(
"../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
"ST000000/SE000001/"
)
# `Path.glob` selects the given relative pattern in the directory represented by the given path
# and yields all matching files (of any type)
# Since `Path.glob` returns a generator, it is converted to a list here for easy display
all_files = list(path_to_head_mri.glob("*"))
mri_data = []
for path in all_files:
# Read the individual DICOM files one by one in the list returned by the path generator
data = pydicom.read_file(path)
mri_data.append(data)
tabulation = Form_Generator()
tabulation.heading_printer(
"Reading all DICOM files that match the specified pattern")
statements = [
"""
path_to_head_mri = Path(
"../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
"ST000000/SE000001/"
)
all_files = list(path_to_head_mri.glob("*"))
mri_data = []
for path in all_files:
data = pydicom.read_file(path)
mri_data.append(data)
"""
]
tabulation.statement_generator(statements)
variables = ["str(path_to_head_mri)", "all_files"]
values = [str(path_to_head_mri), str(reprlib_rules.repr(all_files))]
tabulation.variable_generator(variables, values, 1)
expressions = [
"all_files[0].parts",
"all_files[0].parent",
"all_files[0].name",
"all_files[0].suffixes",
"all_files[0].stem",
"len(all_files)",
"mri_data[0]",
"len(mri_data)",
]
results = [
str(all_files[0].parts),
str(all_files[0].parent),
str(all_files[0].name),
str(all_files[0].suffixes),
str(all_files[0].stem),
str(len(all_files)),
str(reprlib_rules.repr(mri_data[0])),
str(len(mri_data)),
]
tabulation.expression_generator(expressions, results, 1)
Reading all DICOM files that match the specified pattern +---------------------------------------------------------+ | Statement | +---------------------------------------------------------+ | path_to_head_mri = Path( | | "../Datasets/An MRI Dicom Data Set of the Head of a | | Normal Male Human Aged 52/" | | "ST000000/SE000001/" | | ) | | all_files = list(path_to_head_mri.glob("*")) | | | | mri_data = [] | | | | for path in all_files: | | data = pydicom.read_file(path) | | mri_data.append(data) | +---------------------------------------------------------+ +-----------------------+-------------------------------------+ | Variable | Value | +-----------------------+-------------------------------------+ | str(path_to_head_mri) | ../Datasets/An MRI Dicom Data Set | | | of the Head of a Normal Male Human | | | Aged 52/ST000000/SE000001 | | all_files | [PosixPath('../Datasets/An MRI | | | Dicom Data Set of the Head of a | | | Normal Male Human Aged | | | 52/ST000000/SE000001/MR000000'), | | | PosixPath('../Datasets/An MRI | | | Dicom Data Set of the Head of a | | | Normal Male Human Aged | | | 52/ST000000/SE000001/MR000007'), | | | PosixPath('../Datasets/An MRI | | | Dicom Data Set of the Head of a | | | Normal Male Human Aged | | | 52/ST000000/SE000001/MR000009'), | | | PosixPath('../Datasets/An MRI | | | Dicom Data Set of the Head of a | | | Normal Male Human Aged | | | 52/ST000000/SE000001/MR000008'), | | | PosixPath('../Datasets/An MRI | | | Dicom Data Set of the Head of a | | | Normal Male Human Aged | | | 52/ST000000/SE000001/MR000006'), | | | PosixPath('../Datasets/An MRI | | | Dicom Data Set of the Head of a | | | Normal Male Human Aged | | | 52/ST000000/SE000001/MR000001'), | | | ...] | +-----------------------+-------------------------------------+ +-----------------------+-------------------------------------+ | Expression | Result | +-----------------------+-------------------------------------+ | all_files[0].parts | ('..', 'Datasets', 'An MRI Dicom | | | Data Set of the Head of a Normal | | | Male Human Aged 52', 'ST000000', | | | 'SE000001', 'MR000000') | | all_files[0].parent | ../Datasets/An MRI Dicom Data Set | | | of the Head of a Normal Male Human | | | Aged 52/ST000000/SE000001 | | all_files[0].name | MR000000 | | all_files[0].suffixes | [] | | all_files[0].stem | MR000000 | | len(all_files) | 27 | | mri_data[0] | Dataset.file_meta | | | ------------------------------- | | | (0002, 0000) File Meta Information | | | Group Length UL: 214 | | | (0002, 0001) Fil...Group Length | | | UL: 131084 | | | (7fe0, 0010) Pixel Data | | | OW: Array of 131072 | | | elements | | len(mri_data) | 27 | +-----------------------+-------------------------------------+
# DICOM files may not be sorted according to their actual image positions, this can
# be verified by checking Slice Location (0020,1041)
# Slice Location (0020,1041) identifies the relative position of the image plane in millimeters
unordered_slices = {
index: float(slice.SliceLocation) for index, slice in enumerate(mri_data)
}
# For better viewing and processing of 3D data stored with multiple 2D DICOM files,
# these files must be sorted; otherwise the complete scanned image will become cluttered
# and useless
mri_data_ordered = sorted(mri_data, key=lambda slice: slice.SliceLocation)
ordered_slices = {
index: float(slice.SliceLocation) for index, slice in enumerate(mri_data_ordered)
}
tabulation = Form_Generator()
tabulation.heading_printer("Sorting of read DICOM files")
statements = [
"""
unordered_slices = {
index: float(slice.SliceLocation) for index, slice in enumerate(mri_data)
}
mri_data_ordered = sorted(mri_data, key=lambda slice: slice.SliceLocation)
ordered_slices = {
index: float(slice.SliceLocation) for index, slice in enumerate(mri_data_ordered)
}
"""
]
tabulation.statement_generator(statements)
variables = ["unordered_slices", "ordered_slices"]
values = [
str(unordered_slices),
str(ordered_slices),
]
tabulation.variable_generator(variables, values, 1)
expressions = [
"mri_data_ordered[0]",
"len(mri_data_ordered)",
"len(unordered_slices)",
"len(ordered_slices)",
]
results = [
str(reprlib_rules.repr(mri_data_ordered[0])),
str(len(mri_data_ordered)),
str(len(unordered_slices)),
str(len(ordered_slices)),
]
tabulation.expression_generator(expressions, results, 1)
Sorting of read DICOM files +-----------------------------------------------------------+ | Statement | +-----------------------------------------------------------+ | unordered_slices = { | | index: float(slice.SliceLocation) for index, slice in | | enumerate(mri_data) | | } | | | | mri_data_ordered = sorted(mri_data, key=lambda slice: | | slice.SliceLocation) | | | | ordered_slices = { | | index: float(slice.SliceLocation) for index, slice in | | enumerate(mri_data_ordered) | | } | +-----------------------------------------------------------+ +------------------+------------------------------------------+ | Variable | Value | +------------------+------------------------------------------+ | unordered_slices | {0: 0.0, 1: 41.9999963629367, 2: | | | 53.9999958207213, 3: 47.9999970362677, | | | 4: 35.9999959546749, 5: | | | 5.99999663091323, 6: 137.999998321624, | | | 7: 143.999998928727, 8: | | | 71.9999961590453, 9: 89.9999955528687, | | | 10: 83.9999967682912, 11: | | | 77.999996227574, 12: 149.999999502083, | | | 13: 131.999997780749, 14: | | | 23.9999946081714, 15: 17.9999979772582, | | | 16: 11.9999973042441, 17: | | | 29.9999952815023, 18: 107.999995419197, | | | 19: 119.999996566542, 20: | | | 95.9999960937442, 21: 65.9999961939969, | | | 22: 59.9999962290673, 23: | | | 101.999994745866, 24: 125.999997173645, | | | 25: 155.999992554172, 26: | | | 113.999995959439} | | ordered_slices | {0: 0.0, 1: 5.99999663091323, 2: | | | 11.9999973042441, 3: 17.9999979772582, | | | 4: 23.9999946081714, 5: | | | 29.9999952815023, 6: 35.9999959546749, | | | 7: 41.9999963629367, 8: | | | 47.9999970362677, 9: 53.9999958207213, | | | 10: 59.9999962290673, 11: | | | 65.9999961939969, 12: 71.9999961590453, | | | 13: 77.999996227574, 14: | | | 83.9999967682912, 15: 89.9999955528687, | | | 16: 95.9999960937442, 17: | | | 101.999994745866, 18: 107.999995419197, | | | 19: 113.999995959439, 20: | | | 119.999996566542, 21: 125.999997173645, | | | 22: 131.999997780749, 23: | | | 137.999998321624, 24: 143.999998928727, | | | 25: 149.999999502083, 26: | | | 155.999992554172} | +------------------+------------------------------------------+ +-----------------------+-------------------------------------+ | Expression | Result | +-----------------------+-------------------------------------+ | mri_data_ordered[0] | Dataset.file_meta | | | ------------------------------- | | | (0002, 0000) File Meta Information | | | Group Length UL: 214 | | | (0002, 0001) Fil...Group Length | | | UL: 131084 | | | (7fe0, 0010) Pixel Data | | | OW: Array of 131072 | | | elements | | len(mri_data_ordered) | 27 | | len(unordered_slices) | 27 | | len(ordered_slices) | 27 | +-----------------------+-------------------------------------+
# Fill the 3D array in a slice-per-slice manner
full_volume = [slice.pixel_array for slice in mri_data_ordered]
full_volume = np.array(full_volume)
tabulation = Form_Generator()
tabulation.heading_printer(
"One-time extraction of the overall 3D pixel array for all slices"
)
statements = [
"""
full_volume = [slice.pixel_array for slice in mri_data_ordered]
full_volume = np.array(full_volume)
"""
]
tabulation.statement_generator(statements)
variables = ["full_volume"]
values = [str(reprlib_rules.repr(full_volume))]
tabulation.variable_generator(variables, values)
expressions = ["full_volume.shape", "full_volume.dtype"]
results = [str(full_volume.shape), str(full_volume.dtype)]
tabulation.expression_generator(expressions, results)
One-time extraction of the overall 3D pixel array for all slices +-----------------------------------------------+ | Statement | +-----------------------------------------------+ | full_volume = [slice.pixel_array for slice in | | mri_data_ordered] | | full_volume = np.array(full_volume) | +-----------------------------------------------+ +-------------+------------------------------------+ | Variable | Value | +-------------+------------------------------------+ | full_volume | array([[[0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | ..., | | | [0,... ..., | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0]]], | | | dtype=uint16) | +-------------+------------------------------------+ +-------------------+----------------+ | Expression | Result | +-------------------+----------------+ | full_volume.shape | (27, 256, 256) | | full_volume.dtype | uint16 | +-------------------+----------------+
def sort_checker(list_1):
list_2 = list_1[1:]
if all(a <= b for a, b in zip(list_1, list_2)):
return "ascending"
elif all(a >= b for a, b in zip(list_1, list_2)):
return "descending"
elif all(a == b for a, b in zip(list_1, list_2)):
return "identical"
else:
return "unordered"
def title_adder(func):
@wraps(func)
def wrapper(inputs_1, ax, dict_slices, **kwargs):
is_sorted = sort_checker(list(dict_slices.values()))
title = f"DICOM image grid with {is_sorted} slice locations"
func(inputs_1, ax, dict_slices, **kwargs)
ax.set_title(
title,
loc="center",
pad=10,
)
return wrapper
def grid_annotator(func):
@wraps(func)
def wrapper(inputs_1, ax, dict_slices, fontsize_adjustment=0, **kwargs):
x_positions_1, x_positions_2, y_positions, cmap = func(
inputs_1, ax, **kwargs)
for (y, x_1), key in zip(
itertools.product(y_positions, x_positions_1), dict_slices.keys()
):
text = f"No.{(key + 1):>03}"
ax.text(
x_1,
y,
text,
transform=ax.transData,
fontsize=7 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
for (y, x_2), value in zip(
itertools.product(y_positions, x_positions_2), dict_slices.values()
):
text = f"{value:>5.1f} mm"
ax.text(
x_2,
y,
text,
transform=ax.transData,
fontsize=7 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
return wrapper
@title_adder
@grid_annotator
def grid_DICOM_2D_image(inputs_1, ax, n=None, row_size=9, pad_value=1):
if n is None:
n = inputs_1.shape[0]
nrows = math.ceil(n / row_size)
cmap = "bone"
# The only types supported by `torch.Tensor` are float64, float32, float16, complex64,
# complex128, int64, int32, int16, int8, uint8 and bool, not uint16
inputs_1 = torch.Tensor(inputs_1.astype(int)).unsqueeze(1)
row = iter(np.array_split(inputs_1[:n], nrows))
x_positions_1 = [
math.ceil(
x * inputs_1[:row_size].shape[3]
+ 0.075 * inputs_1[:row_size].shape[3]
+ pad_value * x
)
for x in range(inputs_1[:row_size].shape[0])
]
x_positions_2 = [
math.ceil(
x * inputs_1[:row_size].shape[3]
+ 0.55 * inputs_1[:row_size].shape[3]
+ pad_value * x
)
for x in range(inputs_1[:row_size].shape[0])
]
y_positions = [
math.ceil(
y * inputs_1[:row_size].shape[2]
+ 0.125 * inputs_1[:row_size].shape[2]
+ pad_value * y
)
for y in range(nrows)
]
for i in range(nrows):
row_images = next(row)
# Unlike the images in the MNIST dataset, the DICOM file images do not take values
# in the range (0, 1), so the parameter `normalize` needs to be set to True in order to
# move their values proportionally to the range (0, 1)
# By querying the source code of `vision/torchvision/utils.py`, it is clear that any
# single-channel image is converted to a three-channel image, and all are represented
# as tensors, as shown below:
# > if tensor.dim() == 2: # single image H x W
# > tensor = tensor.unsqueeze(0)
# > if tensor.dim() == 3: # single image
# > if tensor.size(0) == 1: # if single-channel, convert to 3-channel
# > tensor = torch.cat((tensor, tensor, tensor), 0)
# > tensor = tensor.unsqueeze(0)
# >
# > if tensor.dim() == 4 and tensor.size(1) == 1: # single-channel images
# > tensor = torch.cat((tensor, tensor, tensor), 1)
grid_row = make_grid(
row_images, nrow=row_size, normalize=True, pad_value=pad_value
)
scalar_row = grid_row.numpy().transpose(1, 2, 0)[:, :, 0]
if i == 0:
scalar_image = scalar_row
else:
scalar_image = np.concatenate((scalar_image, scalar_row), axis=0)
ax.imshow(scalar_image, cmap=cmap)
ax.set(xticks=[], yticks=[], frame_on=False)
ax.grid(False)
return x_positions_1, x_positions_2, y_positions, cmap
unordered_full_volume = [slice.pixel_array for slice in mri_data]
unordered_full_volume = np.array(unordered_full_volume)
descending_slices = {
index: float(slice.SliceLocation)
for index, slice in enumerate(mri_data_ordered[::-1])
}
fig, axs = plt.subplots(3, 1, figsize=(figure_size[0], figure_size[1] / 7 * 9))
grid_DICOM_2D_image(unordered_full_volume, axs[0], unordered_slices)
grid_DICOM_2D_image(full_volume, axs[1], ordered_slices)
grid_DICOM_2D_image(full_volume[::-1], axs[2], descending_slices)
fig.suptitle(
"Visual Comparison of DICOM Image Grids with Unordered and Ordered Slice Locations",
fontsize="x-large",
x=0.5,
y=0,
)
plt.tight_layout()
plt.show()
fig, axs = plt.subplots(9, 3, figsize=(
figure_size[0] / 3 * 2, figure_size[1] / 2 * 7))
for slice_counter, (i, j) in enumerate(itertools.product(range(9), range(3))):
DICOM_2D_image(
axs[i][j],
"bone",
mri_data_ordered[slice_counter],
f"DICOM image with {mri_data_ordered[slice_counter].SliceLocation:^.1f} mm "
"slice location",
fontsize_adjustment=-3,
patien_age="52",
)
fig.suptitle(
"Visualization of MRI Medical Images with Information Displayed in Order of Slice Location",
fontsize="x-large",
x=0.5,
y=0,
)
plt.tight_layout()
plt.show()
# For some image formats such as DICOM, the images also contain associated metadata,
# which is not loaded by the reader by default for time-saving reasons
# `ImageSeriesReader` class provides the ability to read a series of image files into a
# SimpleITK image, and once a series of images has been read, the metadata can be accessed
# directly from the reader
# `GetGDCMSeriesIDs` provides the ability to get all seriesIDs from a DICOM dataset
# The path object must be stringified before using `GetGDCMSeriesIDs`
series_ids = sitk.ImageSeriesReader.GetGDCMSeriesIDs(str(path_to_head_mri))
# `GetGDCMSeriesFileNames` generates a series of file names from the directory containing
# the DICOM dataset and series ID
series_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(
str(path_to_head_mri), series_ids[0]
)
tabulation = Form_Generator()
tabulation.heading_printer("Identification of DICOM image file series")
statements = [
"""
series_ids = sitk.ImageSeriesReader.GetGDCMSeriesIDs(
str(path_to_head_mri)
)
series_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(
str(path_to_head_mri), series_ids[0]
)
"""
]
tabulation.statement_generator(statements)
variables = ["series_ids", "series_file_names"]
values = [
str(reprlib_rules.repr(series_ids)),
str(reprlib_rules.repr(series_file_names)),
]
tabulation.variable_generator(variables, values, 1)
expressions = ["len(series_ids)", "len(series_file_names)"]
results = [str(len(series_ids)), str(len(series_file_names))]
tabulation.expression_generator(expressions, results)
Identification of DICOM image file series +-------------------------------------------------------+ | Statement | +-------------------------------------------------------+ | series_ids = sitk.ImageSeriesReader.GetGDCMSeriesIDs( | | str(path_to_head_mri) | | ) | | series_file_names = | | sitk.ImageSeriesReader.GetGDCMSeriesFileNames( | | str(path_to_head_mri), series_ids[0] | | ) | +-------------------------------------------------------+ +-------------------+---------------------------------------+ | Variable | Value | +-------------------+---------------------------------------+ | series_ids | ('1.3.46.67058...1413262801702',) | | series_file_names | ('../Datasets/...0001/MR000000', | | | '../Datasets/...0001/MR000001', | | | '../Datasets/...0001/MR000002', | | | '../Datasets/...0001/MR000003', | | | '../Datasets/...0001/MR000004', | | | '../Datasets/...0001/MR000005', ...) | +-------------------+---------------------------------------+ +------------------------+--------+ | Expression | Result | +------------------------+--------+ | len(series_ids) | 1 | | len(series_file_names) | 27 | +------------------------+--------+
series_reader = sitk.ImageSeriesReader()
series_reader.SetFileNames(series_file_names)
image_data = series_reader.Execute()
path_to_head_mri
PosixPath('../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/ST000000/SE000001')